\chapter{Electrostatic potentials, forces and torques}
\label{app:force}
In the following, the derivation of forces and torques from the electrostatic potentials employed in our model is presented. %Obviously, the calculation of these quantities is necessary to perform MD simulation. 
General treatments on the derivation of forces and torques from anisotropic potentials can be found elsewhere.~\cite{price84a}%~\cite{cheun76a,price84a,popel94a}

\section{Electrostatic interactions}

\subsection{Charge-charge (Coulombic) interactions}

\subsubsection{Charge-charge potential}
The interaction potential energy $u^{QQ}$ between two point charges $Q_i$ and $Q_j$ located with a distance $r$ between them is defined by Coulomb's law:
\begin{equation}~\label{eq:qq}
u^{QQ}(r)=\frac{Q_iQ_j}{4\pi\epsilon_0r}
\end{equation}
with $\epsilon_0$ the dielectric constant in free space (vacuum).
\subsubsection{Charge-charge forces}
The electrostatic force corresponding to the potential of Equation~\ref{eq:qq} is:
\begin{equation}
\mathbf{f}_{ij}(r)=\frac{Q_iQ_j}{4\pi\epsilon_0r^3}\mathbf{r}
\end{equation}
with $\mathbf{r}=\mathbf{r}_i-\mathbf{r}_j$.

%\subsubsection{Cutoff treatment} To speed up the computation, non-bonded interactions are typically skipped outside a cutoff radius~$r_\mathrm{c}$. In \brahms we have implemented the shifted-force cutoff scheme. 
\paragraph{Shifted-force cutoff scheme} The discontinuity at the cutoff distance $r_\mathrm{c}$ affects both the apparent energy conservation and the actual atomic motion. This discontinuity can be smeared out by changing the form of the potential function slightly, adding a small linear term so that its derivative is zero at the cutoff distance:~\cite{allen,rapa}%\cite[Section~3.3.2] 
\begin{equation}\label{eq:sf} u^{\mathrm{SF}}(r)=u(r)-u_\mathrm{c}-(r-r_\mathrm{c})\frac{\de u(r)}{\de r}\rvert_{r=r_\mathrm{c}}  \end{equation} %The force goes smoothly to zero at the cutoff $r_\mathrm{c}$, removing problems in energy conservation and any numerical instability in the equations of motion.~\cite{allen} %\cite[Section~5.2.4]
For the Coulomb potential, defining a multiplicative shifting function: 
\begin{equation}
S^M(r)=\left(1-\frac{r}{r_c}\right)^2\end{equation}
the Coulomb energy is modified as:
\begin{equation}
u(r)=\frac{q_iq_j}{r}S^M(r)\end{equation}
%There is a consensus that this shifting scheme works well~\cite{levit95a}.
and the force becomes:
\begin{equation}
\mathbf{f}=\frac{q_iq_j}{r^2}\left(1-\frac{r}{r_c}\right)\left(\frac{1}{r}+\frac{1}{r_c}\right)\mathbf{r}=\frac{q_iq_j}{r^3}\left(1-\frac{r}{r_c}\right)\left(1+\frac{r}{r_c}\right)\mathbf{r}=\frac{q_iq_j}{r}\left(\frac{1}{r^2}-\frac{1}{r_c^2}\right)\mathbf{r}
\end{equation}



\subsection{Charge-dipole interaction}

Interaction between a dipole $\mu$ with associated orientation vector $\mathbf{e}_\mu$ and a charge $Q$. In the following we will assume site $j$ to be a point-dipole and site $i$ a point-charge. The angle $\theta$ between the dipole orientation vector  $\mathbf{e}_\mu$  and the interparticle distance vector $\mathbf{r}$ can be obtained from the relation $\cos\theta = \mathbf{e}_\mu\cdot\mathbf{r} / r$. Also, we define $C = Q\mu / 4\pi\epsilon_0$. 

\subsubsection{Charge-dipole potential}
Considering site $j$ a dipole $\mu$, with corresponding orientation vector $\mathbf{e}_\mu$, and site $i$ a charge $Q$, the electrostatic interaction energy is:
\begin{equation}
u^{Q\mu}=\frac{C}{r^3}\,(\mathbf{e}_\mu\cdot\mathbf{r})=\frac{C}{r^2}\,\cos\theta
\end{equation}
with $\mathbf{r}=\mathbf{r}_i-\mathbf{r}_j=\mathbf{r}_Q-\mathbf{r}_\mu$.

\subsubsection{Charge-dipole force}

Pair force:
\begin{equation}
\label{eq:dipChargeForce}
\mathbf{f}_{ij}=-\mathbf{f}_{ji}= \frac{C}{r^3}\left(\frac{3\cos\theta}{r}\mathbf{r}-\hat{\mathbf{e}}\right) = \frac{C}{r^5}\left[3(\hat{\mathbf{e}}_\mu\mathbf{r})\mathbf{r}-r^2\hat{\mathbf{e}}_\mu\right] 
\end{equation}
with $\mathbf{r}=\mathbf{r}_i-\mathbf{r}_j=\mathbf{r}_Q-\mathbf{r}_\mu$.


\subsubsection{Charge-dipole torque}

Pair torque on site $j$:
\begin{equation}
\mathbf{T}_{ji}=\frac{C}{r^3}\left(\mathbf{r}\times\mathbf{e}_{j}\right)=\mathbf{T}_{\mu Q}=\frac{C}{r^3}(\mathbf{r}\times\mathbf{e}_{\mu})
\end{equation}
with $\mathbf{r}=\mathbf{r}_i-\mathbf{r}_j=\mathbf{r}_Q-\mathbf{r}_\mu$.

\subsection{Charge-dipole interaction: shifted-force form}
Considering a charge $Q$ and a dipole $\mu$, being $r=\mathbf{r}_Q-\mathbf{r}_\mu$, the interaction energy is:
\begin{equation}
u^{\mathrm{SF}}(r)=\frac{C}{r^2}\left[1-\left(\frac{r}{r_\mathrm{c}}\right)^2+2\,\frac{r-r_\mathrm{c}}{r_\mathrm{c}}\left(\frac{r}{r_\mathrm{c}}\right)^2\right]\cos\theta=
\frac{C}{r^3}\left[1-3\,\left(\frac{r}{r_\mathrm{c}}\right)^2+2\,\left(\frac{r}{r_\mathrm{c}}\right)^3\right]\mathbf{e}\cdot\mathbf{r}
\end{equation}
where $\cos\theta=\mathbf{e}\cdot\mathbf{r}$.\\
Using the chain rule:
\begin{equation}
\mathbf{f}=-\frac{\partial u}{\partial r}\,\nabla_{\mathbf{r}}r-\frac{\partial u}{\partial \cos\theta}\,\nabla_{\mathbf{r}}\cos\theta
\end{equation}
\begin{equation}
\frac{\partial u}{\partial r}= -\frac{2C\cos\theta}{r^3}\left[1-\left(\frac{r}{r_\mathrm{c}}\right)^3\right]
\end{equation}
Finally we obtain:
\begin{equation}
\mathbf{f}_{\mu Q}=\frac{C}{r^3}\left\{3\,\left[1-\left(\frac{r}{r_\mathrm{c}}\right)^2\right]\frac{\mathbf{r}}{r}\cos\theta-\left[1-3\,\left(\frac{r}{r_\mathrm{c}}\right)^2+2\,\left(\frac{r}{r_\mathrm{c}}\right)^3\right]\mathbf{e}\right\} %= \frac{2C}{r^3}\left[1-\left(\frac{r}{r_\mathrm{c}}\right)^3\right]\mathbf{e}
\end{equation}
\begin{equation}
\mathbf{T}_{\mu Q}=\frac{C}{r^3}\left[1-3\,\left(\frac{r}{r_\mathrm{c}}\right)^2+2\,\left(\frac{r}{r_\mathrm{c}}\right)^3\right](\mathbf{r}\times\mathbf{e})
\end{equation}
Consistency check: for $r_\mathrm{c}\rightarrow\infty$ we recover the unshifted case.

\subsection{Charge-dipole interaction: linear-switch form}
The switched potential energy is:
\begin{equation}
u^{\mathrm{S}}(r)=u(r)\cdot s(r)
\end{equation}
The modulating linear function $s$ is:
\begin{displaymath}
s(r)= 
\left\{ 
\begin{array}{lll}
1 &\mathrm{ if } \qquad r < r_{s}\\
(r_{c}-r)/(r_{c}-r_{s}) &\mathrm{ if }  \qquad r_{s} \leqslant r \leqslant r_{c}\label{eq:ssdmodu} \\
0 &\mathrm{ if } \qquad r > r_{c}\\
\end{array}
 \right. 
\end{displaymath}
The switched force (for $ r_{s} \leqslant r \leqslant r_{c}$) is: 
\begin{equation}
%\mathbf{f}=\frac{C(2r_c-r)}{r^3(r_c-r_s)}\mathbf{e}= \frac{C}{r^3}\left[\left(\frac{3r_\mathrm{c}-2r}{r_\mathrm{c}-r_\mathrm{s}}\right)\frac{\mathbf{r}}{r}\cos\theta-\left(\frac{r_\mathrm{c}-r}{r_\mathrm{c}-r_\mathrm{s}}\right)\mathbf{e}\right]
%\boxed{\mathbf{f}= \frac{C}{r^3}\left[\left(\frac{3r_\mathrm{c}-2r}{r_\mathrm{c}-r_\mathrm{s}}\right)\frac{\mathbf{r}}{r}\cos\theta-\left(\frac{r_\mathrm{c}-r}{r_\mathrm{c}-r_\mathrm{s}}\right)\mathbf{e}\right]}
\mathbf{f}= \frac{C}{r^3(r_\mathrm{c}-r_\mathrm{s})}\left[\left(\frac{3r_\mathrm{c}}{r} - 2\right)\cos\theta\,\mathbf{r}+(r-r_\mathrm{c})\,\mathbf{e}\right]
\end{equation}
The switched torque (for $ r_{s} \leqslant r \leqslant r_{c}$) is:
\begin{equation}
\mathbf{T}=\frac{C(r_c-r)}{r^3(r_c-r_s)}\mathbf{r}\times\mathbf{e}
\end{equation}

\subsection{Dipole-dipole interactions}

Here we consider a pair of (symmetric molecules and, in particular) dipoles $i$ and $j$. The orientations are defined by the unit vectors $\mathbf{e}_i$ and $\mathbf{e}_j$; the angles between these orientation vectors and the interparticle separation vector $\mathbf{r}_{ij}$ are defined respectively as $\theta_i$ and $\theta_j$. Also, we define $\gamma_{ij}$ as the angle between %the plane containing 
$\mathbf{e}_i$ and 
%$\mathbf{r}_{ij}$  and the plane containing 
$\mathbf{e}_j$, %and  $\mathbf{r}_{ij}$, 
so that $\cos{\gamma}_{ij}=\mathbf{e}_i\cdot\mathbf{e}_j$. 
%as in Figure~\ref{fig:orient}.
%\begin{figure} \centering \includegraphics[scale=1]{appendices/orient} \caption[Relative orientation of two dipoles]{The relative orientation of two vectors associated with dipolar sites. From Allen and Tildesley.~\cite{allen}} \label{fig:orient} \end{figure}

\subsubsection{Dipolar potential}
The electrostatic interaction potential energy is:
\begin{equation}
u_{ij}^{\mu\mu}=\frac{\mu^2}{r^3}(\cos\gamma_{ij}-3\cos\theta_i\cos\theta_j)
\end{equation}
where $r = |\mathbf{r}_{ij}|$. Cosines are computed through:
\begin{equation}
\cos\theta_i = \frac{\mathbf{e}_{i}\cdot\mathbf{r}_{ij}}{r|\mathbf{e}_{i}|}\;\; \;\;\;\cos\theta_j = \frac{\mathbf{e}_{j}\cdot\mathbf{r}_{ij}}{r|\mathbf{e}_{j}|}
\end{equation}


\subsubsection{Dipolar forces}

Pair force:
\begin{equation}
\mathbf{f}_{ij}=-\mathbf{f}_{ji}=\frac{3\mu^2}{r^4}\left[(\cos\gamma_{ij}-5\cos\theta_i\cos\theta_j)(\mathbf{r}_{ij}/r+\cos\theta_j\mathbf{e}_i+\cos\theta_i\mathbf{e}_j)\right]
\end{equation}

\subsubsection{Dipolar torques}

Pair torques:
\begin{equation}
\mathbf{T}_{ij}=-\frac{\mu^2}{r^3}\left[\mathbf{e}_i\times\mathbf{e}_j-3\cos\theta_j(\mathbf{e}_i\times\mathbf{r}_{ij})/r\right]
\end{equation}
\begin{equation}
\mathbf{T}_{ji}=-\frac{\mu^2}{r^3}\left[\mathbf{e}_j\times\mathbf{e}_i-3\cos\theta_i(\mathbf{e}_j\times\mathbf{r}_{ij})/r\right]
\end{equation}
\\\\A complete treatment of the dipolar potential, along with the explicit derivation to obtain forces and torques, can be found elsewhere.~\cite{allen} %[p~332-334] 

